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Abstract The macroscopic mechanical properties of 
colloidal particle gels strongly depend on the local ar- 
rangement of the powder particles. Experiments have 
shown that more heterogeneous microstructures exhibit 
up to one order of magnitude higher elastic properties 
than their more homogeneous counterparts at equal vol- 
ume fraction. In this paper, packings of spherical par- 
ticles are used as model structures to computationally 
investigate the elastic properties of coagulated parti- 
cle gels as a function of their degree of heterogeneity. 
The discrete element model comprises a linear elastic 
contact law, particle bonding and damping. The simu- 
lation parameters were calibrated using a homogeneous 
and a heterogeneous microstructure originating from 
earlier Brownian dynamics simulations. A systematic 
study of the elastic properties as a function of the de- 
gree of heterogeneity was performed using two sets of 
microstructures obtained from Brownian dynamics sim- 
ulation and from the void expansion method. Both sets 
cover a broad and to a large extent overlapping range 
of degrees of heterogeneity. The simulations have shown 
that the elastic properties as a function of the degree of 
heterogeneity are independent of the structure genera- 
tion algorithm and that the relation between the shear 
modulus and the degree of heterogeneity can be well 
described by a power law. This suggests the presence 
of a critical degree of heterogeneity and, therefore, a 
phase transition between a phase with finite and one 
with zero elastic properties. 
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1 Introduction 

Random sphere packings are ubiquitous model systems 
for the study of the structural and mechanical proper- 
ties of granular matter or colloids. Geomechanics [I], 
granular flow [2] and mixing and segregation of granu- 
lar materials [3] are but a few examples. The mechan- 
ical properties of granular systems depend on various 
parameters such as the volume fraction [3], the parti- 
cle size distribution [5], material properties as, for ex- 
ample, the particles' friction coefficients [B] or adhesive 
forces [7]. Furthermore, as predicted in [F] and discussed 
in more detail in the following, the mechanical prop- 
erties strongly depend on the microstructure, i.e., the 
local arrangement of the particles. 

The influence of the microstructure on the mechani- 
cal properties is often observed implicitly in experimen- 
tal and computational mechanical tests on structures 
differing in preparation history. Macroscopic stress pro- 
files, for example, where found to strongly depend on 
the sample preparation procedure and thus its mi- 
crostructure [5]. However, systematic investigations of 
the mechanical properties as a function of the mi- 
crostructure are scarce, for two reasons: first, a system- 
atic study of the microstructure-dependent mechanical 
properties requires the possibility of an unambiguous 
characterization of the microstructural arrangement of 
the particles. Second, and equally important, it relies 
on the possibility of a reproducible generation of mi- 
crostructures with distinct local arrangements of the 
particles at constant volume fraction. 
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In [TU], the concept of the degree of heterogene- 
ity (DOH) has been introduced, constituting a scalar 
measure of the heterogeneity of the microstructural ar- 
rangement of a sphere packing. Three distinct struc- 
ture characterization methods in conjunction with pa- 
rameters in fit functions or integrals were shown to al- 
low for a clear quantification and thus classification of 
the DOH of particle structures. In contrast to distri- 
butional structure characterization methods, such as 
the pair correlation function [TT], the bond angle dis- 
tribution the common neighbors distribution |12) 
or the Minkowski functionals [13], the concept of the 
DOH does not rely on the interpretation of distribution 
curves. This quantifiability is seen as an advantage in 
view of a correlation of the mechanical to the structural 
properties. 

Experimentally, the reproducible control of the 
DOH of colloidal microstructures with volume fractions 
between 0.2 and 0.6 is achieved using direct coagulation 
casting [14,15], which is an in situ enzyme-catalyzed 
destabilization method. It allows for the coagulation of 
electrostatically stabilized colloidal suspensions to stiff 
particle structures by either shifting the pH of the sus- 
pension to the particles' isoelectric point or by increas- 
ing the ionic strength of the suspension. Shifting the pH 
leads to "more homogeneous" microstructures through 
diffusion-limited aggregation (ApH-method) while in- 
creasing the ionic strength results in "more heteroge- 
neous" microstructures via reaction-rate-limited aggre- 
gation (Al-method). These differences in heterogeneity 
have been observed using various experimental tech- 
niques such as diffusing wave spectroscopy [16], static 
light transmission 16 or cryogenic scanning electron 
microscopy |17j . 

Rheological and uniaxial compression experiments 
on coagulated colloidal particle structures obtained by 
direct coagulation casting have revealed that those 
with a more heterogeneous microstructure have signif- 
icantly higher elastic moduli than their more homoge- 
neous counterparts [TBI] 191120| . The rheological proper- 
ties, which are the subject of the computational part 
of this study, were investigated experimentally using a 
Bohlin rheometer (Model CS-50, Bohlin Instruments, 
Sweden) equipped with a measuring tool of platc/platc 
geometry (rough surface, 25 mm plate diameter). Os- 
cillatory measurements were performed at a fixed fre- 
quency of 1 Hz with increasing strain amplitude. Fig- 
ure [T] summarizes the measured plateau storage moduli 
G' p , i.e. the elastic shear modulus in the linear regime 
(small deformations), of alumina particle suspensions 
(average particle diameter do — 0.4 /Ltm) destabilized 
by the ApH- and the Al-method, respectively, as a 
function of the volume fraction. For this system, ap- 
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Fig. 1 Experimental plateau storage modulus G' p of alumina 
particle suspensions (average particle diameter do = 0.4 fim) 
formed by the Z\pH- and the Z\I-method of the direct coagulation 
casting process in dependence of the volume fraction I19| . 



proximately four times higher elastic properties are 
measured for heterogeneous than for homogenous mi- 
crostructures at corresponding volume fractions [TW51"] . 

A variation of the ApH-method allowing for a con- 
trolled introduction of heterogeneities is the use of 
alkali-swellable polymer particles. Small amounts of 
these particles, 80 nm in diameter, are admixed to 
the suspended powder particles under acidic condi- 
tions. The polymer particles swell upon changing pH 
during the internal gelling reaction of the direct co- 
agulation casting process and enfold to 0.7 /im in di- 
ameter, thereby rearranging the powder particles and 
thus producing more heterogeneous microstructures. 
These more heterogeneous samples exhibit much higher 
mechanical properties in comparison to samples with- 
out polymer particles. In particular, they present com- 
parably high mechanical properties as samples with 
heterogeneous microstructures produced by the Al- 
method [551155]. 

In summary, strong evidence is given that the dif- 
ferences in macroscopic mechanical properties of coag- 
ulated particle suspensions are controlled by the dif- 
ferences in heterogeneity. A yet unanswered question 
is how these microstructural differences on the length 
scale of a few particle diameters can have such a dra- 
matic influence on the macroscopic mechanical proper- 
ties. 

The aim of this study is to perform a systematic 
computational analysis of the elastic shear properties of 
particle packings at constant volume fraction of <!> = 0.4 
and to correlate these properties with the structures' 
DOH. 
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2 Materials and Methods 

In a first part of this section, the discrete element 
method, which is the simulation method used through- 
out this study, is introduced. Second, the methods, by 
which the two sets of initial microstructures are ob- 
tained are briefly presented. Then, the method and fit 
function that is used in order to quantify the structures' 
DOH is explained and finally, the simulation setup is 
presented. 



2.1 Discrete Element Method 

The discrete element method [23] is a numerical method 
that allows simulating the motion of large numbers of 
interacting particles. The simulations presented in this 
paper are performed using the particle flow code in 
three dimensions PFC 3D [25] . which is an implemen- 
tation of the discrete element method using spherical 
particles. One calculation step consists in the determi- 
nation of the total force on each particle and the nu- 
merical integration of the equations of motion using a 
central difference scheme. The force calculation and the 
determination of the next positions and velocities are 
then repeated until the end of the simulation. 

Our model comprises a linear elastic contact law 
between particles, particle bonding and damping. The 
linear elastic contact law relates the contact forces act- 
ing on two particles in contact linearly to the relative 
displacement between the particles. In PFC 3D a soft- 
contact approach is used, wherein the rigid particles are 
allowed to overlap at contact points. The magnitude of 
the normal contact force F n is given by 

F n = k n U n , (1) 

where k n denotes the normal stiffness and U n is the 
overlap distance (U n > 0). The shear stiffness k s relates 
incremental displacements in shear direction AU S to the 
shear contact force AF S via 

AF S - k s AU s . (2) 

The linear elastic contact law is thus parameterized by 
the normal and shear stiffnesses. 

The so-called contact bonds used in PFC 3D essen- 
tially extend the linear elastic contact law, equations (JJJ 
and @, to particles that are separated, i.e., to parti- 
cles without overlap (U n < 0). These bonds are charac- 
terized by a normal and shear strengths, F„ and Ff . 
A bond breaks if either the normal or the shear bond 
strength is exceeded in normal and shear direction, re- 
spectively. 



Table 1 Simulation parameters. 



Parameter 


Symbol 


Value 


Number of particles 


N 


8 000 


Particle radius 


ro 


0.25 /urn 


Normal particle stiffness 




50 - 150 N/m 


Shear particle stiffness 


Kg 


5-15 N/m 


Normal bond strength 


FS 


10" 4 N 


Shear bond strength 


F S B 


10" 6 N 


Damping coefficient 


a 


0.9 


Friction coefficient 


H 


0.0 


Volume fraction 


<P 


0.4 


Particle density 


P 


3 690 kg/m 3 



Energy dissipation is simulated via a local, non- 
viscous damping term added to the equations of mo- 
tion [26 . This damping force Fd is characterized by 
the damping coefficient a and is given componentwise 
(index i G {1,2,3}) by 

F dti = -a\Fi\siga(vi), (3) 

where F is the force acting on a particle, v is the par- 
ticle's velocity and 

f +1, if x > 
sign(x) = < -1, if x < . (4) 
[ 0, ifx = Q 

The damping force is thus proportional to the force 
acting on the particle and indeed leads to energy- 
dissipation since Fd • v < 0. In particular, only acceler- 
ating motion is damped |25j . An analogous expression 
to Eq. ([3]) is used to describe the damping of the parti- 
cles' rotational motion. 

In conclusion, the simulation model is characterized 
by five microscopic parameters: the stiffness of the par- 
ticles and the bond strength both for normal and shear 
direction, and the damping coefficient. These parame- 
ters, in the following referred to as microparameters, are 
compiled in Table [TJ In particular, the normal to shear 
particle stiffness ratio k n /k s = 10 is kept constant and a 
value of 0.9 is chosen for the damping coefficient, which 
corresponds to the values used in previous studies [TOj 
|2"T] . In PFC 3D the frictional term between two parti- 
cles is superseded in the presence of a bond. However, 
it may become active in the case of larger deformation 
when unbonded particles come into contact. In order 
to prevent such contacts to result in friction, a friction 
coefficient /i = was used. It has to be mentioned that 
the present paper does not include a detailed sensitivity 
study of every simulation parameter. Here, the focus is 
on the influence of the microstructural differences on 
the mechanical properties. This influence can a priori 
also be observed with other choices of parameters as 
long as a unique set is used for all structures. 
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2.2 Initial Microstructures 

Two sets of microstructures are used as initial particle 
configurations for the simulation of the elastic prop- 
erties. A first set originates from previous Brownian 
dynamics (BD) simulations [TTll2"8] where the coagula- 
tion of electrostatically stabilized colloidal suspensions 
to stiff particle structures was simulated. The presence 
and depth of a secondary minimum in the inter-particle 
potential, described by the Derjaguin-Landau-Verwey- 
Overbeek theory [25], was shown to account for the 
variations in the DOH of the resulting particle struc- 
tures and is adjusted via the surface potential <?o- For 
<Pb = m V, the electrostatic double layer repulsion is 
zero and the inter-particle potential is only given by 
the attractive van der Waals potential. For the set of 
simulation parameters used in 11111281 , a secondary min- 
imum appears for ^Pq > 12 mV and an energy barrier 
between the local maximum and the secondary mini- 
mum emerges. For \Pq = 15 mV, a repulsive barrier of 
5.65 ksT is present. The model further contains the fric- 
tional Stokes' drag force and a random Brownian force 
caused by the suspending liquid. In this paper, BD- 
microstructures generated with <Po = 0, 12, 13, 14 and 
15 mV are used, exhibiting an increasing DOH with 
increasing tf'o [10] - 

In particular, the most and least heterogeneous BD- 
microstructures have been shown to nicely correspond 
to experimental silica structures using the pair correla- 
tion function [T7]. These two structures are thus used 
here to calibrate the microparameters in order for the 
simulations to reproduce the experimental values given 
in Fig.H](cf. Sec. EE]). 

The second set of initial microstructures is obtained 
using the void expansion method (VEM) presented 
in [3D]. In contrast to BD, which simulates the phys- 
ical processes during coagulation according to estab- 
lished laws and methods, VEM is a purely stochas- 
tic method inspired by the experimental generation 
of heterogeneous microstructures using alkali-swellable 
polymer particles. This method employs so-called void- 
particles that mimic the polymer particles. The void- 
particles, having a small initial diameter, are ran- 
domly placed into the simulation box containing the 
structure-particles. The core part of VEM is the re- 
peated increase in the void-particles' diameter. During 
this procedure, the structure-particles are rearranged 
and pushed into contact. Finally, the void- particles are 
deleted. As shown in [27], VEM allows generating mi- 
crostructures presenting a broad range of DOH, which 
is controlled by the void- to structure-particle num- 
ber ratio. In particular, VEM-microstructures cover 
an approximately 20% larger range of DOH, slightly 



shifted to higher values than the BD-microstructures. 
In this study, the number of void-particles Ny ranges 
between 1 000 and 16 000. The void-particles' normal 
and shear stiffness used in VEM are k n y = 10 2 N/m 
and k s y = 10~ 2 N/m, respectively. 

All structures are submitted to the same relaxation 
procedure before undergoing the simulated mechanical 
testing. Bonds between neighboring particles having a 
ccnter-to-center separation smaller than d e = (1 + e)2ro 
with e = 0.01 are installed. Then, several calculation 
steps are performed, where, after each step, the trans- 
lational and rotational velocities are set to zero. This 
allows for an efficient reduction of particle overlaps and 
thus of the internal strain energy in the structures with- 
out significant changes in the particle positions. The re- 
laxation process is aborted as soon as the mean unbal- 
anced force in the structure has decreased to a constant 
value. 

The structures investigated in this study have equal 
volume fraction of <P = 0.4, consist of N = 8 000 
monodispersed spherical particles having a radius ro = 
0.25 /xm and are contained in a cubic simulation box 
with periodic boundaries and edge length Lb ox given 

by 



J box 



''o 



/4Vtt\ 



1/3 



(5) 



The structures used in this study are provided as elec- 
tronic supplementary material (Online Resources 1-12). 

2.3 Degree of Heterogeneity 

In [TU] , three distinct methods allowing for a quantifica- 
tion of the DOH using scalar measures were introduced. 
Here, the DOH is characterized by means of the cumu- 
lative pore size distribution P(rp > r) estimated using 
the exclusion probability [3T]. Equation © was shown 
to nicely fit P(r P > r) of both the VEM- and the BD- 
microstructures [T0ll2"7] . 



P(r P > r) = 1 - erf 



r/r - b 
aV2 



(6) 



The DOH is quantified by the width of the error 
function measured by parameter a. The values of a for 
the various microstructures are summarized in Table [2] 
Parameter b represents the location of the maximum 
of the underlying Gaussian distribution, i.e., the most 
probable pore to particle radius ratio. 

Additionally, the structures' coordination numbers 
C are given in Tableland plotted in Fig. [2] against the 
DOH-parameter a. This figure shows that the DOH and 
the coordination number are to a large extent indepen- 
dent. 
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Fig. 2 Degree of heterogeneity as measured by parameter a as 
a function of the coordination number C for all structures used 
in this study. 
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Oscillation Cycle 

Fig. 4 Stress crjg (open squares, left scale) and strain 7jg (full 
circles, right scale) as a function of shear oscillation cycle. 
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Fig. 3 Rheological simulation setup. 

2.4 Simulation Setup 

The simulation setup is schematically shown in Fig. [21 
The structure (white particles) is placed between peri- 
odic boundaries in x- and y-direction. In z-direction, the 
boundary conditions arc imposed via three additional 
particle layers (gray), obtained from the periodic rep- 
etition of the structure. The lower layers, representing 
the lower plate in the shear experiment, are immobile 
whereas the upper layers shear the sample periodically 
along the x-axis with a fixed frequency of 1 Hz and 
constant volume. The shear amplitude is increased af- 
ter each oscillation. In particular, ten oscillations with a 
linearly increasing deflection of the upper plate up to a 
maximum deflection of 0.1% Lb ox are performed. Spher- 
ical regions, so-called measurement spheres [25) . are de- 
fined in different heights monitoring the mean stress 



tensor in their respective region. The absolute value 
of the complex modulus G is obtained by dividing the 
peak shear stress value with the actual peak strain of 
each oscillation. The storage modulus G' is then calcu- 
lated by multiplication with the cosine of the phase an- 
gle between excitation and response, which corresponds 
to the evaluation of G' according to 32 . In total, the 
simulations are performed six times for each structure. 
The shear plate's normal vector is successively oriented 
along the x-, y- and z-axis and, for each orientation, the 
structure is sheared along the two remaining axes. 



3 Results and Discussion 

3.1 Calibration of the Model 

Typical stress and strain curves are shown in Fig. HJ 
presenting the shear stress <j\i (open squares, left scale) 
and the applied shear deflection 712 (full circles, right 
scale) as a function of the shear oscillation cycle. o\2 
increases linearly with 712 and no phase shift between 
excitation and response is observed (phase angle < 
10~ 4 7r). This confirms that the simulations present a 
purely elastic behavior and the constant ratio 012/712 
is identified with the plateau storage modulus G' p . 

The particle stiffness k n has been calibrated using 
the most homogeneous and the most heterogeneous BD- 
microstructures (i'o = mV and Wo = 15 mV, respec- 
tively) using a constant ratio k n /k s = 10. Its influence 
on the mechanical properties of these two structures is 
shown in Fig. [5] (left graph). The experimental values 
are shown on the right. 

The simulated values for G' p quantitatively agree 
with the experimental values within the experimental 
error for normal stiffness values k n ranging between 50 
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Fig. 5 Simulated plateau storage moduli G' p of the most and 
least heterogeneous BD-microstructure as a function of the par- 
ticle normal stiffness k n . Right: Experimental G^-values for alu- 
mina particle structures (cf. Fig. [![). ^ = OA for all microstruc- 
turcs. 

and 150 N/m. For both microstructures, a linear depen- 
dency is found between G' p and k n , which is expressed in 
a constant ratio between the G p -values for the heteroge- 
neous and homogeneous microstructure as a function of 
k n : G' p HE /G' p HO fn 3.4. This is in good agreement with 
the average experimental G^-ratio of 3.6 at <!> = 0.4. 
For all further simulations the particle normal stiffness 
is fixed to k n = 100 N/m. 

The inter-particle bond strength was strictly speak- 
ing not part of the calibration process. In order to fo- 
cus on the initial elastic properties, its value was cho- 
sen high enough to prevent any bond breakage for the 
deformations considered in this study. The chosen nor- 
mal bond strength F„ = 10~ 4 N, in conjunction with 
the particle normal stiffness k n — 100 N/m allows for a 
maximum particle separation distance of 10 -6 m, which 
corresponds to two particle diameters. Given the maxi- 
mum shear displacement of 0.001 L\, OXl where L\, ox is of 
the order of roughly 20 particle diameters, this value is 
never exceeded. 



3.2 Elastic Properties as a Function of the DOH 

Using a particle normal stiffness of 100 N/m, the 
plateau storage moduli of the various VEM- and BD- 
microstructure were simulated. The resulting G' - values 
are summarized in Table [5] and are shown in Fig. [5] as 
a function of the structures' DOH expressed by param- 
eter a. Triangles and circles correspond to BD- and 
VEM-microstructures, respectively. Each data point 
in Fig. [5] corresponds to one structure sheared in all 
six directions. The error bars are thus reflecting the 
anisotropy of the structures. 



Table 2 DOH-parameter a I10II27| . coordination number C and 
simulated plateau storage moduli G' p of the various VEM- and 
BD-microstructures. 





Structure 


a 


C 


G' p in MPa 




Ny 


= 1000 


0.594 


4.68 


1.73 ±0.383 




N v 


= 2 000 


0.521 


4.71 


1.83 ±0.136 




N v 


= 4 000 


0.467 


4.69 


1.62 ±0.190 


VEM 


N v 


= 8 000 


0.418 


4.69 


1.32 ±0.237 




N v 


= f 3 000 


0.399 


4.69 


1.32 ±0.267 




N v 


= f 6 000 


0.392 


4.61 


1.17 ±0.231 




&a -- 


= mV 


0.373 


4.52 


0.480 ± 0.049 




V -- 


= 12 mV 


0.404 


4.84 


1.19 ±0.141 


BD 


V -- 


= 13 mV 


0.410 


4.75 


1.05 ±0.124 




<P -- 


= 14 mV 


0.472 


5.08 


1.46 ±0.120 




W -- 


= 15 mV 


0.536 


4.97 


1.62 ±0.213 



2.2 i- 
2.0 - 




0.35 a g 0.40 0.45 0.50 0.55 0.60 

a 



Fig. 6 Simulated plateau storage moduli G' p of the various mi- 
crostructures (<P = 0.4, rg = 0.25 fiva) as a function of the DOH 
in terms of parameter a and power law fit (solid line). 

Two conclusions can be drawn from Fig. [S] first, 
the elastic properties present a clear dependence on the 
microstructure's heterogeneity such that G' p increases 
for increasing DOH-parameter a. Second, the behav- 
ior is independent of the structure generation algo- 
rithm. Indeed, the elastic moduli of the VEM- and the 
BD-microstructures agree within the error bars for mi- 
crostructures with comparable DOH. 

The solid line in Fig. [6] represents a power law fit 
given by 

G' p oc {a-a Q f,a > a , (7) 

with a,Q — 0.373 ±0.001 (shown as dashed line in Fig.[B]) 
and = 0.187 ± 0.032. The use of this fit function 
is inspired by percolation theory where a power law 
scaling is found for the elastic properties as a function 
of the volume fraction [331I34] . In view of this analogy 
and based on the simulation data currently available, 
this fit function suggests a phase transition with ao the 
critical DOH. Below this value, the elastic properties 
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are zero; above this value, they increase for increas- 
ing a. A DOH-value a = 0.265 below ao was found for 
the stabilized colloidal microstructures with repulsive 
inter-particle potentials [lOj . In the model used in this 
work, this microstructure would indeed present negli- 
gible elastic properties since there would be no bonds 
between the particles. 

It has to be mentioned that percolating structures 
presenting a DOH below the critical value exist. Such 
structures, as for example intermediate steps during the 
simulated coagulation toward the most homogeneous 
microstructure using BD, are not expected to have zero 
elastic properties. However, they have an average coor- 
dination number significantly below the range consid- 
ered in this study and are therefore not included. Also, 
it has to be emphasized that the obtained value for 
do does not solely rely on the data point possessing the 
lowest a-value. Indeed, a power law fit without this data 
point yields a' — 0.360, which is roughly 4% below do- 

The uncertainty of the value for the critical DOH ao 
roots in fitting the power law (Eq. to the simulated 
plateau storage moduli of the samples studied here us- 
ing error propagation. The limits for ao do therefore 
not represent an uncertainty obtained by analyzing a 
statistically significant control sample. A more exten- 
sive study, involving a larger number of microstructures 
possessing a DOH in the immediate vicinity of ao, i.e. 
between 0.35 and 0.4, would be required in order to 
determine ao more accurately. 

4 Summary and Conclusions 

In this paper, the relation between the microstruc- 
ture and the elastic properties of packings of spher- 
ical particles was investigated computationally using 
the discrete element method. Two distinct sets of ini- 
tial microstructures were subjected to strain-controlled 
oscillatory shear simulations. The first set originated 
from Brownian dynamics simulations [Tl"ll2"8] where the 
physical processes during coagulation were simulated. 
The second set was obtained using the void expansion 
method |30j . In contrast to Brownian dynamics simula- 
tions, the latter is a purely stochastic method. Both sets 
cover a broad range of degrees of heterogeneity, which, 
in this work, was characterized using the width of the 
pore size distribution. 

The simulations of the mechanical properties have 
been performed using a discrete element model with five 
microparameters: the particle stiffness and the bond 
strength in normal and shear direction, respectively 
and damping. The inter-particle normal and shear bond 
strengths were set to a high value in order to pre- 
vent any bond breaking for the deformations applied 



during the shear simulations. This allowed focussing 
on the initial, purely elastic behavior found in exper- 
iment [20]. The particle stiffness was calibrated using 
the two Brownian dynamics microstructures presenting 
the highest and lowest DOH, which correspond to ex- 
perimental silica structures [17] . 

A quantitative agreement with experiment was 
achieved. For particle normal stiffnesses ranging be- 
tween 50 and 150 N/m, the absolute values of the 
experimental elastic moduli are reproduced, as shown 
in Fig. [5] In particular, the ratio between the elastic 
moduli of the heterogeneous and the homogeneous mi- 
crostructures is in good agreement with the experimen- 
tal value. This latter result is particularly remarkable, 
since the G^-ratio has not been subject of any calibra- 
tion but reflects the influence of the microstructures and 
emphasizes the strong influence of the DOH indepen- 
dently of the interaction potential. Indeed, the model 
used in this study constitutes an important simplifi- 
cation with respect to theoretical models such as the 
Derjaguin-Landau-Verwey-Overbeek theory [21] or the 
Johnson-Kendall- Roberts theory [35] . The former gives 
the inter-particle potential as the sum of the van der 
Waals attraction and the electrostatic repulsion and the 
latter describes the adhesive force between particles. 

Using the calibrated model, the elastic properties 
of all microstructures generated using the void expan- 
sion method and Brownian dynamics were simulated. 
This investigation has revealed a correlation between 
the elastic properties and the structures' DOH. It has 
also shown that the elastic properties are independent 
of the structure generation algorithm. Indeed, for struc- 
tures with comparable degrees of heterogeneity, com- 
parable plateau storage moduli were found. A power 
law fit was shown to well reproduce the relation be- 
tween the plateau storage modulus and the DOH above 
a critical value, which indicates a phase transition be- 
tween a phase with finite and a phase with zero elastic 
properties. This result is in nice analogy to earlier find- 
ings 33 , 34 , where the elastic properties were shown to 
exhibit a power law dependence as a function of the vol- 
ume fraction. This study thus suggests an extension to 
percolation theory, where usually the volume fraction 
constitutes the continuous variable presenting a critical 
value. Here, the DOH, possessing the same critical be- 
havior as the volume fraction, has been identified as an 
additional variable. 

The results presented in this paper show that the 
degree of heterogeneity is particularly useful in order to 
quantify and characterize the heterogeneity of particle 
packings since it allows for an unprecedented theoret- 
ical description of their elastic properties. However, it 
does not directly lead to a better understanding of the 
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mechanisms, by which the elastic properties increase for 
increasing heterogeneity. Such an understanding must, 
in some way, include an intermediate step, such as an 
analysis of the distribution of particle contacts using 
the fabric tensor, for example (36) . or of load-bearing 
substructures, as discussed in [T2] . 
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